!  IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in znaupd .
!  MODE, ISHIFT, MXITER: see the definition of IPARAM in znaupd .
!
!  NP      Integer.  (INPUT/OUTPUT)
!          Contains the number of implicit shifts to apply during
!          each Arnoldi iteration.
!          If ISHIFT=1, NP is adjusted dynamically at each iteration
!          to accelerate convergence and prevent stagnation.
!          This is also roughly equal to the number of matrix-vector
!          products (involving the operator OP) per Arnoldi iteration.
!          The logic for adjusting is contained within the current
!          subroutine.
!          If ISHIFT=0, NP is the number of shifts the user needs
!          to provide via reverse comunication. 0 < NP < NCV-NEV.
!          NP may be less than NCV-NEV since a leading block of the current
!          upper Hessenberg matrix has split off and contains "unwanted"
!          Ritz values.
!          Upon termination of the IRA iteration, NP contains the number
!          of "converged" wanted Ritz values.
!
!  IUPD    Integer.  (INPUT)
!          IUPD == 0: use explicit restart instead implicit update.
!          IUPD /= 0: use implicit update.
!
!  V       complex(8)  N by (NEV+NP) array.  (INPUT/OUTPUT)
!          The Arnoldi basis vectors are returned in the first NEV
!          columns of V.
!
!  LDV     Integer.  (INPUT)
!          Leading dimension of V exactly as declared in the calling
!          program.
!
!  H       complex(8)  (NEV+NP) by (NEV+NP) array.  (OUTPUT)
!          H is used to store the generated upper Hessenberg matrix
!
!  LDH     Integer.  (INPUT)
!          Leading dimension of H exactly as declared in the calling
!          program.
!
!  RITZ    complex(8)  array of length NEV+NP.  (OUTPUT)
!          RITZ(1:NEV)  contains the computed Ritz values of OP.
!
!  BOUNDS  complex(8)  array of length NEV+NP.  (OUTPUT)
!          BOUNDS(1:NEV) contain the error bounds corresponding to
!          the computed Ritz values.
!
!  Q       complex(8)  (NEV+NP) by (NEV+NP) array.  (WORKSPACE)
!          Private (replicated) work array used to accumulate the
!          rotation in the shift application step.
!
!  LDQ     Integer.  (INPUT)
!          Leading dimension of Q exactly as declared in the calling
!          program.
!
!  WORKL   complex(8)  work array of length at least
!          (NEV+NP)**2 + 3*(NEV+NP).  (WORKSPACE)
!          Private (replicated) array on each PE or array allocated on
!          the front end.  It is used in shifts calculation, shifts
!          application and convergence checking.
!
!
!  IPNTR   Integer array of length 3.  (OUTPUT)
!          Pointer to mark the starting locations in the WORKD for
!          vectors used by the Arnoldi iteration.
!          -------------------------------------------------------------
!          IPNTR(1): pointer to the current operand vector X.
!          IPNTR(2): pointer to the current result vector Y.
!          IPNTR(3): pointer to the vector B * X when used in the
!                    shift-and-invert mode.  X is the current operand.
!          -------------------------------------------------------------
!
!  WORKD   complex(8)  work array of length 3*N.  (WORKSPACE)
!          Distributed array to be used in the basic Arnoldi iteration
!          for reverse communication.  The user should not use WORKD
!          as temporary workspace during the iteration !!!!!!!!!!
!          See Data Distribution Note in ZNAUPD .
!
!  RWORK   real(8)    work array of length  NEV+NP ( WORKSPACE)
!          Private (replicated) array on each PE or array allocated on
!          the front end.
!
!  INFO    Integer.  (INPUT/OUTPUT)
!          If INFO == 0, a randomly initial residual vector is used.
!          If INFO /= 0, RESID contains the initial residual vector,
!                          possibly from a previous run.
!          Error flag on output.
!          =     0: Normal return.
!          =     1: Maximum number of iterations taken.
!                   All possible eigenvalues of OP has been found.
!                   NP returns the number of converged Ritz values.
!          =     2: No shifts could be applied.
!          =    -8: Error return from LAPACK eigenvalue calculation;
!                   This should never happen.
!          =    -9: Starting vector is zero.
!          = -9999: Could not build an Arnoldi factorization.
!                   Size that was built in returned in NP.
!-----------------------------------------------------------------------
subroutine znaup2(ido, which, nev, np, tol, resid, ishift, mxiter, v, n, h, ncv, workl, ipntr, workd, info)
    implicit none
    character  which*2
    integer    ido, info, ishift, ncv, n, mxiter, nev, np, ipntr(3)
    real(8)    tol
    complex(8) h(ncv, ncv), resid(n), v(n, ncv), workd(3 * n), workl(ncv * (ncv + 3))
    complex(8), parameter :: zero = (0.0D0, 0.0D0)
    ! Local Scalars
    complex(8), external :: zdotc
    real(8), external :: dznrm2, dlamch, dlapy2
    logical cnorm, getv0, initv, update, ushift
    integer ierr, iter, nconv, nevbef, nev0, np0, nptemp, kplusp, i, j
    complex(8) ritz(ncv), bounds(ncv), q(ncv, ncv)
    real(8) rnorm, eps23, rtemp
    character  wprime*2
    save       cnorm, getv0, initv, update, ushift, rnorm, iter, kplusp, nconv, nevbef, nev0, np0, eps23

    if (ido == 0) then
        nev0 = nev
        np0 = np
        ! kplusp is the bound on the largest Lanczos factorization built. nconv is the current number
        ! of "converged" eigenvalues. iter is the counter on the current iteration step.
        kplusp = nev + np
        nconv = 0
        iter = 0
        ! Get machine dependent constant.
        eps23 = dlamch('Epsilon-Machine')
        eps23 = eps23**(2.0D0 / 3.0D0)
        ! Set flags for computing the first NEV steps of the Arnoldi factorization.
        getv0 = .true.
        update = .false.
        ushift = .false.
        cnorm = .false.
        if (info /= 0) then
            ! User provides the initial residual vector.
            initv = .true.
            info = 0
        else
            initv = .false.
        end if
    end if
    ! Get a possibly random starting vector and force it into the range of the operator OP.
    if (getv0) then
        call zgetv0(ido, initv, n, 1, v, resid, rnorm, workd, info)
        if (ido /= 99) return

        if (rnorm == 0.0D0) then
            ! The initial vector is zero. Error exit.
            info = -9
            mxiter = iter
            nev = nconv
            ido = 99
            return
        end if
        getv0 = .false.
        ido = 0
    end if
    ! Back from reverse communication : continue with update step
    if (update) go to 20
    ! Back from computing user specified shifts
    if (ushift) go to 50
    ! Back from computing residual norm at the end of the current iteration
    if (cnorm)  go to 100
    ! Compute the first NEV steps of the Arnoldi factorization
    call znaitr(ido, n, 0, nev, resid, rnorm, v, h, ncv, ipntr, workd, info)
    if (ido /= 99) return
    if (info > 0) then
        np = info
        mxiter = iter
        info = -9999
        ido = 99
        return
    end if

    ! MAIN  ARNOLDI  ITERATION  LOOP
    ! Each iteration implicitly restarts the Arnoldi factorization in place.
    do while (.TRUE.)
        iter = iter + 1
        ! Compute NP additional steps of the Arnoldi factorization.
        ! Adjust NP since NEV might have been updated by last call
        ! to the shift application routine znapps.
        np = kplusp - nev
        ! Compute NP additional steps of the Arnoldi factorization.
        ido = 0
20      continue
        update = .true.
        call znaitr(ido, n, nev, nev + np, resid, rnorm, v, h, ncv, ipntr, workd, info)

        if (ido /= 99) return

        if (info > 0) then
            np = info
            mxiter = iter
            info = -9999
            ido = 99
            return
        end if
        update = .false.

        ! Compute the eigenvalues and corresponding error bounds of the current upper Hessenberg matrix.
        call zneigh(rnorm, kplusp, h, ritz, bounds, q, ierr)
        if (ierr /= 0) then
            info = -8
            ido = 99
            return
        end if
        ! Select the wanted Ritz values and their bounds to be used in the convergence test.
        ! The wanted part of the spectrum and corresponding error bounds are in the last NEV loc.
        ! of RITZ, and BOUNDS respectively.
        nev = nev0
        np = np0
        ! Make a copy of Ritz values and the corresponding Ritz estimates obtained from zneigh.
        call zcopy(kplusp, ritz, 1, workl(kplusp**2 + 1), 1)
        call zcopy(kplusp, bounds, 1, workl(kplusp**2 + kplusp + 1), 1)
        ! Select the wanted Ritz values and their bounds to be used in the convergence test.
        ! The wanted part of the spectrum and corresponding
        ! bounds are in the last NEV loc. of RITZ BOUNDS respectively.
        call zsortc(which, nev + np, ritz, bounds)
        if (ishift == 1) then
            call zsortc('SM', np, bounds, ritz)
        end if
        ! Convergence test: currently we use the following criteria.
        ! The relative accuracy of a Ritz value is considered acceptable if:
        ! error_bounds(i) <= tol*max(eps23, magnitude_of_ritz(i)).
        nconv = 0
        do i = 1, nev
            rtemp = max(eps23, dlapy2(dble(ritz(np + i)), dimag(ritz(np + i))))
            if (dlapy2(dble(bounds(np + i)), dimag(bounds(np + i))) <= tol * rtemp) then
                nconv = nconv + 1
            end if
        end do
        ! Count the number of unwanted Ritz values that have zero |
        ! Ritz estimates. If any Ritz estimates are equal to zero |
        ! then a leading block of H of order equal to at least    |
        ! the number of Ritz values with zero Ritz estimates has  |
        ! split off. None of these Ritz values may be removed by  |
        ! shifting. Decrease NP the number of shifts to apply. If |
        ! no shifts may be applied, then prepare to exit
        nptemp = np
        do j = 1, nptemp
            if (bounds(j) == zero) then
                np = np - 1
                nev = nev + 1
            end if
        end do

        if (nconv >= nev0 .or. iter > mxiter .or. np == 0) then
            ! Prepare to exit. Put the converged Ritz values and corresponding bounds in RITZ(1:NCONV)
            ! and BOUNDS(1:NCONV) respectively. Then sort. Be careful when NCONV > NP
            ! Use h( 3,1 ) as storage to communicate rnorm to zneupd  if needed
            h(3, 1) = dcmplx(rnorm, 0.0D0)
            ! Sort Ritz values so that converged Ritz values appear within the first NEV locations
            ! of ritz and bounds, and the most desired one appears at the front.
            if (which == 'LM') wprime = 'SM'
            if (which == 'SM') wprime = 'LM'
            if (which == 'LR') wprime = 'SR'
            if (which == 'SR') wprime = 'LR'
            if (which == 'LI') wprime = 'SI'
            if (which == 'SI') wprime = 'LI'
            call zsortc(wprime, kplusp, ritz, bounds)
            ! Scale the Ritz estimate of each Ritz value by 1 / max(eps23, magnitude of the Ritz value).
            do j = 1, nev0
                rtemp = max(eps23, dlapy2(dble(ritz(j)), dimag(ritz(j))))
                bounds(j) = bounds(j) / rtemp
            end do
            ! Sort the Ritz values according to the scaled Ritz estimates.  This will push all the converged ones
            ! towards the front of ritz, bounds (in the case when NCONV < NEV.)
            wprime = 'LM'
            call zsortc(wprime, nev0, bounds, ritz)
            ! Scale the Ritz estimate back to its original value.
            do j = 1, nev0
                rtemp = max(eps23, dlapy2(dble(ritz(j)), dimag(ritz(j))))
                bounds(j) = bounds(j) * rtemp
            end do
            ! Sort the converged Ritz values again so that  the "threshold" value appears at the front of ritz and bound.
            call zsortc(which, nconv, ritz, bounds)
            ! Max iterations have been exceeded.
            if (iter > mxiter .and. nconv < nev0) info = 1
            ! No shifts to apply.
            if (np == 0 .and. nconv < nev0)  info = 2
            np = nconv
            mxiter = iter
            nev = nconv
            ido = 99
            return
        else if (nconv < nev0 .and. ishift == 1) then
            ! Do not have all the requested eigenvalues yet. To prevent possible stagnation, adjust the size of NEV.
            nevbef = nev
            nev = nev + min(nconv, np / 2)
            if (nev == 1 .and. kplusp >= 6) then
                nev = kplusp / 2
            else if (nev == 1 .and. kplusp > 3) then
                nev = 2
            end if
            np = kplusp - nev
            ! If the size of NEV was just increased resort the eigenvalues.
            if (nevbef < nev) then
                call zsortc(which, kplusp, ritz, bounds)
                call zsortc('SM', np, bounds, ritz)
            end if
        end if

        if (ishift == 0) then
            ! User specified shifts: pop back out to get the shifts and return them in the first 2*NP locations of WORKL.
            ushift = .true.
            ido = 3
            return
        end if
50      continue
        ushift = .false.
        if (ishift /= 1) then
            ! Move the NP shifts from WORKL to RITZ, to free up WORKL for non-exact shift case.
            call zcopy(np, workl, 1, ritz, 1)
        end if
        ! Apply the NP implicit shifts by QR bulge chasing.
        ! Each shift is applied to the whole upper Hessenberg matrix H.
        ! The first 2*N locations of WORKD are used as workspace.
        call znapps(n, nev, np, ritz, v, n, h, ncv, resid, q, workd)
        ! Compute the B-norm of the updated residual.
        ! Keep B*RESID in WORKD(1:N) to be used in the first step of the next call to znaitr.
        cnorm = .true.
        call zcopy(n, resid, 1, workd, 1)
100     continue
        ! Back from reverse communication; WORKD(1:N) := B*RESID
        rnorm = dznrm2(n, resid, 1)
        cnorm = .false.
    end do
    ! End of main iteration loop
    mxiter = iter
    nev = nconv
    ido = 99
end
